CLL Cytokine Screen 2021: Figure 5
9 November 2021
1 Figure 5
In this sub-vignette we present the analysis and source code for figure 5. This sub-vignette can be built along with all other sub-vignettes, by running CLLCytokineScreen2021.Rmd.
1.1 Set up
Load Libraries
library(patchwork)
library(ggplot2)
library(data.table)
library(magrittr)
library(ggplot2)
library(pheatmap)
library(tidyr)
library(dplyr)Set plot directory
plotDir = ifelse(exists(".standalone"), "", "../../inst/figs/")
if(plotDir!="") if(!file.exists(plotDir)) dir.create(plotDir)1.2 Load data
Raw
#df: tibble containing all screening viability data
load( "../../data/df.RData")#from tsvs
df <- read.table(file= "../../inst/extdata/df.txt",header = TRUE) %>% as_tibble()
df$Cytokine <- factor(df$Cytokine, levels= c("No Cytokine",
"Resiquimod",
"IL-4",
"TGF-b1",
"IL-1b",
"Interferon gamma",
"SDF-1a",
"sCD40L" ,
"sCD40L+IL-4",
"soluble anti-IgM",
"CpG ODN",
"IL-6",
"IL-10",
"IL-21",
"HS-5 CM",
"IL-15",
"BAFF",
"IL-2" ))
df$Drug <- factor(df$Drug, levels =c("DMSO",
"BAY-11-7085",
"Everolimus",
"Fludarabine",
"I-BET 762",
"Ibrutinib","Idelalisib",
"Luminespib",
"Nutlin-3a",
"PRT062607",
"Pyridone-6",
"Ralimetinib",
"Selumetinib"))Process data
# List of Cytokines, Drugs and treatments
##Drugs
thedrugs <- unique(df$Drug)
##Cytokines
thecytokines <- unique(df$Cytokine)
#Treatments
#drugs
drugtreatments <- list()
drugtreatments <-
lapply(thedrugs, function(x){
paste("treatment_drug", x, sep='')
})
#cytokines
cytokinetreatments <- list()
cytokinetreatments <-
lapply(thecytokines, function(y){
paste("treatment_cytokine", y, sep='')
})
# Tidy up data frame to use in linear model
df <-
dplyr::filter(df,
#use high concentrations only
Drug_Concentration %in% c("High","None")) %>%
#select columns for linear model
dplyr::select(PatientID, Drug, Cytokine, Log) %>%
#rename columns
plyr::rename(c("Drug"="treatment_drug", "Cytokine"="treatment_cytokine", "Log"="Viability"))1.3 Set aesthetics
source("../../R/themes_colors.R")1.4 Pre - work
Fit the linear model
# define the base level per treatment_type as the "no"-treatment
df$treatment_drug <- as.factor(df$treatment_drug)
df$treatment_cytokine <- as.factor(df$treatment_cytokine)
df$treatment_drug %<>% relevel("DMSO")
df$treatment_cytokine %<>% relevel("No Cytokine")
# inspect design matrix of interaction model
# model.matrix(~treatment_drug + treatment_cytokine + treatment_drug:treatment_cytokine, df)
# fit model with interaction
fit <- lm(Viability ~ treatment_drug * treatment_cytokine, df)Extract p values
#extract p values
pvalues <- summary(fit)$coefficients[,4]
#order as a dataframe
pvaldf <- as.data.frame(as.matrix(pvalues)) %>%
setDT(keep.rownames = TRUE)
#rename columns
colnames(pvaldf) <- c("treatment", "pvalue")
#filter out single agent treatments
singletreatments <- c(drugtreatments, cytokinetreatments, "(Intercept)")
pvaldf <- dplyr::filter(pvaldf, !treatment %in% singletreatments)
pvaldf <-
pvaldf %>%
#remove "treatment" prefix
mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))
#split drug:cytokine into two columns
pvaldf <-
data.frame(pvaldf, do.call(rbind, strsplit(pvaldf$treatment, split = ":", fixed = TRUE)))
#renove treatment column
pvaldf <- pvaldf[,c("pvalue","X1", "X2")]
#rename columns
colnames(pvaldf) <- c("pvalue","drug", "Cytokine")1.5 Plot Figures
1.5.1 Fig 5A
dummydata <-
as.data.frame(list(c(1,2,3,4),
c(0,0,0,0),
c(-0.6,-0.6,0,-0.1),
c(0.6,0.6,0.1,-0.1),
c(0.5,-0.5,0.8,-0.8)),
col.names=c("plot",
"coord1",
"coord2",
"coord3",
"coord4")) %>%
as_tibble()%>%
mutate(expected = coord2 + coord3, measured = coord4) %>%
pivot_longer(cols = coord1:measured, names_to = c("plot"), values_to="y.value", names_repair="unique" ) %>%
dplyr::rename(Plot = plot...1, x.value = plot...2) %>%
mutate(x.value = factor(x.value, levels = c("coord1", "coord2", "coord3", "coord4", "expected", "measured"))) %>%
#add annotation of interaction type
mutate(Nature=if_else(Plot%in%c(1,2),"Antagonistic", "Synergistic"),
Nature=factor(Nature, levels = c("Antagonistic", "Synergistic")),
Direction=if_else(Plot%in%c(1,3), "Positive \u03B2 int", "Negative \u03B2 int"),
Direction=factor(Direction, levels = c("Positive \u03B2 int", "Negative \u03B2 int")),
Annotation=case_when(Plot==1~"I",
Plot==2~"II",
Plot==3~"III",
Plot==4~"IV"),
Plot=as.character(Plot)
)## New names:
## * plot -> plot...1
## * plot -> plot...2Fig5A <-
ggplot(dplyr::filter(dummydata, !x.value%in%c("expected", "measured")),
aes(x=x.value, y=y.value, group=Plot, color=Plot))+
geom_hline(yintercept = 0, linetype=2, color="black")+
geom_point(size=1)+
geom_line(size=1)+
geom_text(aes(label=Annotation), x=1, y=.8, size=7)+
scale_color_manual(values = c("#A6093D","#003DA5","#A6093D","#003DA5"))+
guides(color="none")+
geom_segment(data=dplyr::filter(dummydata, grepl('expected', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="blue", size=2) +
geom_segment(data=dplyr::filter(dummydata, grepl('measured', x.value)), aes(x = 3.5, y = y.value, xend = 4.5, yend = y.value), color="black", size=2)+
theme_bw() +
t1 +
xlab("") +
ylab("Log(Viability)") +
coord_cartesian(ylim = c(-1, 1)) +
scale_y_continuous(n.breaks=3)+
scale_x_discrete(labels=c( "coord1"="DMSO","coord2"="Drug","coord3"="Stimulus","coord4"="Drug +\nStimulus")) +
theme(axis.ticks.x= element_blank(),axis.ticks.y = element_blank(),
plot.tag=element_text(size = 20)) +
facet_grid(Nature~Direction)
Fig5A
1.5.2 Fig 5B
Barplot quantifying interaction types
fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Condition = rn, comb_coefficient = V1) %>%
dplyr::filter(grepl("drug", Condition)&grepl("cytokine", Condition)) %>%
mutate(Condition = gsub("treatment_drug", "", Condition)) %>%
mutate(Condition = gsub("treatment_cytokine", "", Condition)) %>%
separate(Condition, c("Drug", "Cytokine"), sep=":")
drug_fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Drug = rn, drug_coefficent=V1) %>%
dplyr::filter(grepl("drug", Drug)&!grepl("cytokine", Drug)) %>%
mutate(Drug = gsub("treatment_drug", "", Drug))
cyt_fit_as_df <-
fit$coefficients %>%
as.matrix() %>%
as.data.frame() %>%
setDT(keep.rownames = TRUE) %>%
as_tibble() %>%
dplyr::rename(Cytokine=rn, cyt_coefficent=V1) %>%
dplyr::filter(!grepl("drug", Cytokine)&grepl("cytokine", Cytokine)) %>%
mutate(Cytokine = gsub("treatment_cytokine", "", Cytokine))
fit_combinations <-
left_join(fit_as_df, drug_fit_as_df, by = "Drug") %>%
left_join(cyt_fit_as_df, by = "Cytokine") %>%
mutate(predicted_coefficient=cyt_coefficent+drug_coefficent) %>%
mutate(obs_coef=predicted_coefficient+comb_coefficient) %>%
left_join(mutate(pvaldf,drug=as.character(drug),Cytokine=as.character(Cytokine)), by=c("Drug"="drug", "Cytokine"="Cytokine")) %>%
dplyr::filter(pvalue<=0.05) %>%
mutate(Category_Symbol=case_when(
comb_coefficient>=0&(obs_coef<cyt_coefficent|obs_coef<drug_coefficent)~ "I",
comb_coefficient>=0&(obs_coef>cyt_coefficent&obs_coef>drug_coefficent)~ "III",
comb_coefficient<0&(obs_coef>cyt_coefficent|obs_coef>drug_coefficent)~ "II",
comb_coefficient<0&(obs_coef<cyt_coefficent&obs_coef<drug_coefficent)~ "IV"
)) %>%
mutate(Category=case_when(
Category_Symbol=="I"~ "Positive Antagonistic",
Category_Symbol=="III"~ "Positive Synergistic",
Category_Symbol=="II"~ "Negative Antagonistic",
Category_Symbol=="IV"~ "Negative Synergistic"
)) %>%
mutate( Category=factor(Category, levels=c("Positive Antagonistic","Negative Antagonistic","Positive Synergistic","Negative Synergistic")))
Fig5B <-
fit_combinations %>%
dplyr::group_by(Category) %>%
tally() %>%
ggplot(aes(x=Category, y=n, fill=Category))+
geom_bar(width = 1, stat = "identity")+
t1+
guides(fill="none")+
xlab("")+ ylab("")+
scale_fill_manual(values = c("#A6093D","#003DA5","#A6093D","#003DA5"))
Fig5B
1.5.3 Fig 5C
Heatmap of significant drug : stimulus interaction coefficients
Get interaction coefficients matrix for plotting heatmap
#extract beta interaction values
betavalues <- summary(fit)$coefficients[,1]
#create matrix of beta values
betadf <- as.data.frame(as.matrix(betavalues)) %>%
setDT(keep.rownames = TRUE)
#rename columns
colnames(betadf) <- c("treatment", "betavalue")
#remove beta values for drugs and cytokines alone
betadf <- dplyr::filter(betadf, !treatment %in% singletreatments)
#remove treatment prefix
betadf <- betadf %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_drug", "", .)))) %>%
mutate_at(vars(treatment), funs(as.character(gsub("treatment_cytokine", "", .))))
#split drug:cytokine by colon
betadf <- data.frame(betadf, do.call(rbind, strsplit(betadf$treatment, split = ":", fixed = TRUE)))
#select columns of interest
betadf <- betadf[,c("betavalue","X1", "X2")]
#rename columns
colnames(betadf) <- c("betavalue","drug", "Cytokine")#set significance
a <- 0.05
#make a matrix of beta values, where beta = 0 if corresponding p val is > than significance threshold
pvalmat <- xtabs(pvalue~drug+Cytokine, data=pvaldf)
bvalmat <- xtabs(betavalue~drug+Cytokine, data=betadf)
#check in same order
# bvalmat[,colnames(pvalmat)]
# bvalmat[rownames(pvalmat),]
bvalmat[pvalmat >= a] = 0
#transform matrix
bvalmat <- t(bvalmat)
tree <-
pheatmap(bvalmat, silent = TRUE)
cyt_order <- rownames(bvalmat[tree$tree_row[["order"]],])
drug_order <- colnames(bvalmat[,tree$tree_col[["order"]]])
Fig5C <-
#append p values
left_join(betadf, pvaldf, by = c("drug", "Cytokine")) %>%
#make beta value 0 if p value is < sig
mutate(betavalue = ifelse(pvalue<a, betavalue, 0)) %>%
#put drugs in order of clustered heatmap
mutate(drug=factor(drug, levels = drug_order)) %>%
#put stimuli in order of clustered heatmap
mutate(Cytokine=factor(Cytokine, levels = rev(cyt_order))) %>%
#make heatmap with ggplot
ggplot(aes(x=drug, y=Cytokine))+
geom_tile(aes(fill=betavalue),color = "grey")+
#set colour scale
scale_fill_gradientn(colors=c("#003DA5", "white", "#A6093D"), limits=c(-0.6,0.6))+
#set theme
t2+
theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.ticks = element_blank(),
axis.ticks.x =element_blank(),
panel.background = element_blank(),
axis.line.x = element_blank(),
axis.line.y = element_blank(),
legend.key.height=unit(2.5, "cm"))+
labs(fill = expression(paste(beta["int"], "-value")))+
#add category
geom_text(data = fit_combinations, aes(x=Drug, label=Category_Symbol) )+
scale_y_discrete(labels=c("TGF-b1"="TGF-\u03B21", "sCD40L+IL-4"="sCD40L + IL4", "IL-1b"="IL1\u03B2", "IL-4"="IL4", "IL-6"="IL6","IL-15"="IL15","IL-10"="IL10", "IL-21"="IL21","IL-2"="IL2", "Interferon gamma"= "Interferon \u03B3", "SDF-1a"="SDF-1\u03B1"))
Fig5C
1.5.4 Fig 5D - G
Preparation: Line plots of examples of drug and stimulus combinations that show an interaction
#get lists of drugs and cytokines without baseline treatments
thedrugs.only <- thedrugs %>% setdiff("DMSO") %>% sort()
thecytokines.only <- thecytokines %>% setdiff("No Cytokine")%>% sort()
# get a list of interactions for which p value of beta interaction is significant
sign_conditions <-
pvaldf %>%
dplyr::filter(pvalue<=0.05) %>%
mutate(condition = paste0("treatment_drug:",
drug,
":treatment_cytokine:",
Cytokine))
gg <- vector(mode="list", length=length(sign_conditions$condition))
names(gg) <- sign_conditions$condition
plotList <-
lapply(names(gg), function(i){
#preparations for the interaction plot
drug <- dplyr::filter(sign_conditions, condition ==i)$drug
drugTreat <- paste("treatment_drug", drug, sep='')
cyt <- dplyr::filter(sign_conditions, condition ==i)$Cytokine
cytokineTreat <-paste("treatment_cytokine", cyt, sep='')
Category <- fit_combinations[which(fit_combinations$Drug==drug &
fit_combinations$Cytokine==cyt),]$Category
#extract coefficients to plot
## baseline effect
baseline <- fit$coefficients["(Intercept)"]
## single treatment effect for type_1 = "a"
single_1 <- fit$coefficients[drugTreat]
## single treatment effect for type_2 = "b"
single_2 <- fit$coefficients[cytokineTreat]
# interaction Viability
interaction <- fit$coefficients[paste(drugTreat, cytokineTreat, sep=':')]
## observed single treatment values
observed_single_1 <- baseline + single_1
observed_single_2 <- baseline + single_2
## observed double treatment value
observed_double <- baseline + single_1 + single_2 + interaction
## predicted double treatment value (based only on single effects)
predicted_double <- baseline + single_1 + single_2
xlabels <- c("DMSO", drug,
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
cyt=="soluble anti-IgM"~"soluble\nanti-IgM",
TRUE~cyt),
paste(drug, " +\n",
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
TRUE~cyt),
sep=""))
segment_size <- 2
plotTab <-
df %>%
dplyr::filter(treatment_drug %in% c("DMSO" , drug)) %>%
dplyr::filter(treatment_cytokine %in% c("No Cytokine", cyt)) %>%
dplyr::mutate(treatment_combination =
#if baseline treatment
ifelse(treatment_drug=="DMSO" &
treatment_cytokine=="No Cytokine",
"DMSO",
#if drug only treatment
ifelse(treatment_drug == drug &
treatment_cytokine =="No Cytokine",
paste0("Drug=", drug),
#if cytokine only treatment
ifelse(treatment_drug == "DMSO" &
treatment_cytokine == cyt,
paste0("Cytokine=", cyt),
#otherwise combinatorial treatment
paste0("Drug=", drug, "\nCytokine=", cyt))))) %>%
#make treatment_combination as a factor
mutate(treatment_combination = factor(treatment_combination,
levels=c("DMSO",
paste0("Drug=", drug),
paste0("Cytokine=", cyt),
paste0("Drug=", drug, "\nCytokine=", cyt)
)))
#make plot
ggplot(plotTab, aes(treatment_combination, Viability, group=PatientID)) +
geom_hline(yintercept = 0, linetype=2, color="black")+
lemon::geom_pointline( size=1,
colour=ifelse(interaction > 0,
"#A6093D","#003DA5"), alpha=0.3) +
#add predicted viability with control treatment
geom_segment( size=segment_size,
aes(x=1-.2, xend=1+.2, y=baseline, yend=baseline),
colour = "black") +
#add predicted viability with drug only treatment
geom_segment( size=segment_size,
aes(x=2-.2, xend=2+.2, y=observed_single_1,yend=observed_single_1),
colour="black") +
#add predicted viability with cytokine only treatment
geom_segment( size=segment_size,
aes(x=3-.2, xend=3+.2, y=observed_single_2, yend=observed_single_2),
colour="black") +
#add predicted viability with both treatments
geom_segment( size=segment_size,
aes(x=4-.2, xend=4+.2, y=observed_double, yend=observed_double),
colour="black") +
#add predicted viability with both treatments, without interaction
geom_segment( size=segment_size,
aes(x=4-.2, xend=4+.2, y=predicted_double,
yend=predicted_double), color="blue") +
#add category and combination as title
ggtitle(paste( Category,": ", drug, " & ",
case_when(cyt=="IL-4"~"IL4",
cyt=="Interferon gamma"~"Interferon \u03B3",
cyt=="sCD40L+IL-4"~"sCD40L + IL4",
TRUE~cyt), sep="")) +
#add theme
t2 + theme(axis.text.x = element_text(size=fontsize+5)) +
theme(plot.tag=element_text(size = 30))+
scale_x_discrete(labels = xlabels) +
ylab("Log (Viability)") +
xlab("")
} )
names(plotList) <- names(gg)1.5.4.1 Fig 5D
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)
1.5.4.2 Fig 5E
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)
1.5.4.3 Fig 5F
wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)
1.5.4.4 Fig 5G
wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)
1.5.4.5 Fig 5H
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)
1.5.4.6 Fig 5G
wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)
1.6 Arrange plots
design1 <-"
1223
4455
6677
8899
"
Fig5<-
wrap_elements(Fig5A) + theme(plot.tag=element_text(size = 30)) +
wrap_elements(Fig5B) + theme(plot.tag=element_text(size = 30)) +
wrap_elements(Fig5C) + theme(plot.tag=element_text(size = 30)) +
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:IL-4`)+
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:Interferon gamma`)+
wrap_plots(plotList$`treatment_drug:Pyridone-6:treatment_cytokine:sCD40L+IL-4`)+
wrap_plots(plotList$`treatment_drug:Ralimetinib:treatment_cytokine:Interferon gamma`)+
wrap_plots(plotList$`treatment_drug:Ibrutinib:treatment_cytokine:CpG ODN`)+
wrap_plots(plotList$`treatment_drug:Luminespib:treatment_cytokine:soluble anti-IgM`)+
plot_layout(design=design1,
heights = c(1,0.5,0.5,0.5),
widths = c(0.92,0.08,0.05,0.95)) +
plot_annotation( tag_levels = "A", title="Figure 5", theme = theme(title=element_text(size = 20)))
Fig5
1.7 Appendix
Sys.info()## sysname
## "Darwin"
## release
## "19.6.0"
## version
## "Darwin Kernel Version 19.6.0: Thu May 6 00:48:39 PDT 2021; root:xnu-6153.141.33~1/RELEASE_X86_64"
## nodename
## "mac-huber20.local"
## machine
## "x86_64"
## login
## "holly"
## user
## "holly"
## effective_user
## "holly"sessionInfo()## R version 4.1.0 (2021-05-18)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_1.0.7 tidyr_1.1.3 pheatmap_1.0.12 magrittr_2.0.1
## [5] data.table_1.14.0 ggplot2_3.3.5 patchwork_1.1.1 BiocStyle_2.20.2
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.7 highr_0.9 plyr_1.8.6
## [4] RColorBrewer_1.1-2 pillar_1.6.2 compiler_4.1.0
## [7] BiocManager_1.30.16 tools_4.1.0 digest_0.6.27
## [10] lattice_0.20-44 evaluate_0.14 lifecycle_1.0.0
## [13] tibble_3.1.3 gtable_0.3.0 pkgconfig_2.0.3
## [16] rlang_0.4.11 DBI_1.1.1 magick_2.7.2
## [19] yaml_2.2.1 lemon_0.4.5 xfun_0.25
## [22] gridExtra_2.3 withr_2.4.2 stringr_1.4.0
## [25] knitr_1.33 generics_0.1.0 vctrs_0.3.8
## [28] grid_4.1.0 tidyselect_1.1.1 glue_1.4.2
## [31] R6_2.5.0 fansi_0.5.0 rmarkdown_2.10
## [34] bookdown_0.23 farver_2.1.0 purrr_0.3.4
## [37] scales_1.1.1 ellipsis_0.3.2 htmltools_0.5.1.1
## [40] assertthat_0.2.1 colorspace_2.0-2 labeling_0.4.2
## [43] utf8_1.2.2 stringi_1.7.3 munsell_0.5.0
## [46] crayon_1.4.1